function results = VARmodel(DATA,nlag,c_case,DATA_EX,nlag_ex)
% =======================================================================
% Performs vector autogressive estimation
% =======================================================================
% results = VARmodel(DATA,nlag,c_case,DATA_EX)
% -----------------------------------------------------------------------
% INPUTS 
%   DATA = an (nobs x neqs) matrix of y-vectors
%	nlag = the lag length
%
% OPTIONAL INPUTS
%   c_case  = 0 no constant; 1 with constant ; 2 constant + time trend (default = 1)
%	DATA_EX = optional matrix of variables (nobs x nvar_ex)
%
% OUTPUT
%    results.X         = independent variable (ordered as in VARmakexy)
%    results.Y         = dependent variable
%    results.X_EX      = matrix of exogenous variables
%    results.beta      = matrix of estimated coefficients (ordered as in VARmakexy)
%    results.sigma     = VCV matrix of residuals (adjusted for dof)
%    results.nobs      = nobs, # of observations adjusted (i.e., removing the lags)
%    results.neqs      = neqs, # of equations
%    results.nlag      = nlag, # of lags
%    results.nvar      = nlag*neqs, # of endogenous variables per equation
%    results.nvar_ex   = nvar_ex, # of exogenous variables;
%    results.ntotcoeff = nlag*neqs+nvar_ex+c_case , # coefficients to be estimated per per equation
%--------------------------------------------------------------------------
%    results(eq).beta  = bhat for equation eq
%    results(eq).tstat = t-statistics 
%    results(eq).tprob = t-probabilities
%    results(eq).resid = residuals 
%    results(eq).yhat  = predicted values 
%    results(eq).y     = actual values 
%    results(eq).sige  = e'e/(n-nvar)
%    results(eq).rsqr  = r-squared
%    results(eq).rbar  = r-squared adjusted
%    results(eq).boxq  = Box Q-statistics
%    results(eq).ftest = Granger F-tests
%    results(eq).fprob = Granger marginal probabilities
%
% Note: compared to Eviews, there is a difference in the estimation of the 
% constant when lag is > 2. This is because Eviews initialize the trend
% with the number of lags (i.e., when lag=2, the trend is [2 3 ...T]), 
% while VARmakexy.m initialize the trend always with 1.
% =======================================================================
% Ambrogio Cesa Bianchi, May 2012
% ambrogio.cesabianchi@gmail.com



%% Check inputs
%==============
[nobs, neqs] = size(DATA);

% Check if ther are constant, trend, both, or none
if ~exist('c_case','var')
    c_case = 1;
end

% Check if there are exogenous variables///////////////////////////////////
if exist('DATA_EX','var')
    [nobs2, num_ex] = size(DATA_EX);
    % Check that DATA and DATA_EX are conformable
    if (nobs2 ~= nobs)
        error('var: nobs in DATA_EX-matrix not the same as y-matrix');
    end
    clear nobs2
else
    num_ex = 0;
end

% Check if there is lag orderd of DATA_EX, otherwise set it to 0
if ~exist('nlag_ex','var')
    nlag_ex = 0;
end


%% Save some parameters and create data for VAR estimation
%=========================================================
nobse                 = nobs - max(nlag,nlag_ex);
    results.nobs      = nobse;
    results.neqs      = neqs;
    results.nlag      = nlag;
    results.nlag_ex   = nlag_ex;
nvar                  = neqs*nlag; 
    results.nvar      = nvar;
nvar_ex               = num_ex*(nlag_ex+1);
    results.nvar_ex   = nvar_ex;
ntotcoeff             = nvar + nvar_ex + c_case;
    results.ntotcoeff = ntotcoeff;
    results.c_case    = c_case;

% Create independent vector and lagged dependent matrix
[Y, X] = VARmakexy(DATA,nlag,c_case);

% Create (lagged) exogenous matrix
if nvar_ex
    X_EX  = VARmakelags(DATA_EX,nlag_ex);
    if nlag == nlag_ex
        X = [X X_EX];
    elseif nlag > nlag_ex
        diff = nlag - nlag_ex;
        X_EX = X_EX(diff+1:end,:);
        X = [X X_EX];
    elseif nlag < nlag_ex
        diff = nlag_ex - nlag;
        Y = Y(diff+1:end,:);
        X = [X(diff+1:end,:) X_EX];
    end
end


%% OLS estimation equation by equation
%=====================================

% pull out each y-vector and run regressions
for j=1:neqs
    Yvec = Y(:,j);
    ols_struct = ols(Yvec,X);
    aux = ['eq' num2str(j)];
    eval( ['results.' aux '.beta  = ols_struct.beta;'] );        % bhats
    eval( ['results.' aux '.tstat = ols_struct.tstat;'] );       % t-stats
    % compute t-probs
    tstat = zeros(nvar,1);
    tstat = ols_struct.tstat;
    tout = tdis_prb(tstat,nobse-nvar);
    eval( ['results.' aux '.tprob = tout;'] );                   % t-probs
    eval( ['results.' aux '.resid = ols_struct.resid;'] );       % resids 
    eval( ['results.' aux '.yhat = ols_struct.yhat;'] );         % yhats
    eval( ['results.' aux '.y    = Yvec;'] );                    % actual y
    eval( ['results.' aux '.rsqr = ols_struct.rsqr;'] );         % r-squared
    eval( ['results.' aux '.rbar = ols_struct.rbar;'] );         % r-adjusted
    eval( ['results.' aux '.sige = ols_struct.sige;'] );
    
    % do the Q-statistics
    % use residuals to do Box-Pierce Q-stats
    % use lags = nlag in the VAR
    % NOTE: a rule of thumb is to use (1/6)*nobs but this seems excessive to me
    elag = mlag(ols_struct.resid,nlag);
    % feed the lags
    etrunc = elag(nlag+1:nobse,:);
    rtrunc = ols_struct.resid(nlag+1:nobse,1);
    qres   = ols(rtrunc,etrunc);
    if nlag ~= 1
    	boxq = (qres.rsqr/(nlag-1))/((1-qres.rsqr)/(nobse-nlag));
    else
        boxq = (qres.rsqr/(nlag))/((1-qres.rsqr)/(nobse-nlag));
    end

    eval( ['results.' aux '.boxq = boxq;'] );

    % form matrices for joint F-tests (exclude each variable sequentially)
    for r=1:neqs
        xtmp = [];
        for s=1:neqs
            if s ~= r
                xlag = mlag(DATA(:,s),nlag);
                if nlag == nlag_ex
                    xtmp = [xtmp trimr(xlag,nlag,0)];
                elseif nlag > nlag_ex
                    xtmp = [xtmp trimr(xlag,nlag,0)];
                elseif nlag < nlag_ex
                    xtmp = [xtmp trimr(xlag,nlag+diff,0)];
                end
            end
        end
        % we have an xtmp matrix that excludes 1 variable
        % add deterministic variables (if any) and constant term
        if nvar_ex > 0
            [xtmp X_EX ones(nobse,1)];
        else
            xtmp = [xtmp ones(nobse,1)];
        end
        % get ols residual vector
        b = xtmp\Yvec; % using Cholesky solution
        etmp = Yvec-xtmp*b;
        sigr = etmp'*etmp;
        % joint F-test for variables r
        sigu = ols_struct.resid'*ols_struct.resid;
        ftest(r,1) = ((sigr - sigu)/nlag)/(sigu/(nobse-nvar)); 
    end

    eval( ['results.' aux '.sige = ols_struct.sige;'] );
%     eval( ['results.' aux '.ftest = ftest;' ]);     
%     eval( ['results.' aux '.fprob = fdis_prb(ftest,nlag,nobse-nvar);' ]); 

end % end of loop over equations 

%% Compute the matrix of coefficients & VCV
%==========================================
BETA = (X'*X)\(X'*Y); % ordered by block lags (i.e., all coefficients with first lag, then second,...)
results.beta = BETA;

% %------- ADDED by Tamoni in 2021 -----------------------------------------%
% % To impose exogeneity on the first variable, if needed
% yHHL{1,1}     =  {DATA(:,1)'};
% for ii=1:size(   DATA(:,1:end),2)-1
%     y2{ii,1} =   DATA(:,ii+1)';
% end
% yHHL{2,1}     =  y2;
% LL{1,1}     = nlag; 
% LL{2,1}     = repmat(nlag,size(DATA(:,2:end),2),1);
% trend{1,1}  = 0;
% trend{2,1}  = [zeros(       size(DATA(:,2:end),2),1)];
% [B, S, et] = MLEVARGC(yHHL, LL, trend); 
% results.beta(1:end-num_ex,1) = B(1,[1 3:end])';
% if num_ex>0
% results.beta(end-nlag+1:end,1) = 0;
% end
% %-------------------------------------------------------------------------%

SIGMA = (1/(nobse))*(Y-X*BETA)'*(Y-X*BETA); %(nobse-ntotcoeff) adjusted for # of estimated coeff per equation
results.sigma = SIGMA;
results.residuals = Y - X*BETA;

results.X = X;
results.Y = Y;
if nvar_ex > 0
    results.X_EX = X_EX;
end

if exist('DATA_EX','var')
results.nlagsNL=size(DATA_EX,2);
end







%%

function [Y, X] = VARmakexy(DATA,lags,c_case)
% =======================================================================
% Builds the VAR process from the data-matrix DATA. It orders the data into
% the Y and X matrix --> Example: [x y] = [x(-1) y(-1) x(-2) y(-2)]
% =======================================================================
% [Y, X] = VARmakexy(DATA, lags, c_case)
% -----------------------------------------------------------------------
% INPUT
%   DATA   : matrix containing the original data
%   lags   : lag order of the VAR
%   c_case : 0, no constant, no trend
%            1, constant, no trend
%            2, constant, trend
% OUPUT
%   Y: dependent variable
%   X: independent variable
% =======================================================================
% Ambrogio Cesa Bianchi, July 2011
% ambrogio.cesabianchi@gmail.com


[nobs, m]= size(DATA);

%Y matrix 
Y = DATA(lags+1:end,:);

%X-matrix 
if c_case==0
        X=[];
        for jj=0:lags-1
            X = [DATA(jj+1:nobs-lags+jj,:), X];
        end
        
elseif c_case==1 %constant
        X = [];
        for jj=0:lags-1
            X = [DATA(jj+1:nobs-lags+jj,:), X];
        end
       X = [ones(nobs-lags,1) X];
       
elseif c_case==2 % time trend and constant
        X = [];
        for jj=0:lags-1
            X = [DATA(jj+1:nobs-lags+jj,:), X];
        end
        trend=1:size(X,1);
        X = [ones(nobs-lags,1) trend' X];
     
% elseif c_case==3 % squared time trend, linear time trend, and constant
%         X = [];
%         for jj=0:lags-1
%             X = [DATA(jj+1:nobs-lags+jj,:), X];
%         end
%         trend=1:size(X,1);
%         X = [ones(nobs-lags,1) trend' (trend').^2 X];
end

function X = VARmakelags(DATA,lags)
% =======================================================================
% Builds a matrix with lagged values of DATA, i.e. if DATA = [x y],
% VARmakelags(DATA,1) yields X = [x y x(-1) y(-1)]
% =======================================================================
% X = VARmakelags(DATA,lags)
% -----------------------------------------------------------------------
% INPUT
%   DATA : matrix containing the original data
%   lags : lag order
%
% OUPUT
%   X    : matrix of lagged values
% =======================================================================
% Ambrogio Cesa Bianchi, February 2012
% ambrogio.cesabianchi@gmail.com


[T]= size(DATA,1);

% Create the lagged matrix
X = [];
for jj=0:lags-1
    X = [DATA(jj+1:T-lags+jj,:), X];
end

% Finally, save the non-lagged values...
aux = DATA(lags+1:end,:);

%... and append to the lagged matrix
X = [aux X];


function results=ols(y,x)
% PURPOSE: least-squares regression 
%---------------------------------------------------
% USAGE: results = ols(y,x)
% where: y = dependent variable vector    (nobs x 1)
%        x = independent variables matrix (nobs x nvar)
%---------------------------------------------------
% RETURNS: a structure
%        results.meth  = 'ols'
%        results.beta  = bhat     (nvar x 1)
%        results.tstat = t-stats  (nvar x 1)
%        results.bstd  = std deviations for bhat (nvar x 1)
%        results.yhat  = yhat     (nobs x 1)
%        results.resid = residuals (nobs x 1)
%        results.sige  = e'*e/(n-k)   scalar
%        results.rsqr  = rsquared     scalar
%        results.rbar  = rbar-squared scalar
%        results.dw    = Durbin-Watson Statistic
%        results.nobs  = nobs
%        results.nvar  = nvars
%        results.y     = y data vector (nobs x 1)
%        results.bint  = (nvar x2 ) vector with 95% confidence intervals on beta
%---------------------------------------------------
% SEE ALSO: prt(results), plt(results)
%---------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jlesage@spatial-econometrics.com
%
% Barry Dillon (CICG Equity)
% added the 95% confidence intervals on bhat

if (nargin ~= 2); error('Wrong # of arguments to ols'); 
else
 [nobs nvar] = size(x); [nobs2 junk] = size(y);
 if (nobs ~= nobs2); error('x and y must have same # obs in ols'); 
 end;
end;

results.meth = 'ols';
results.y = y;
results.nobs = nobs;
results.nvar = nvar;

if nobs < 10000
  [q r] = qr(x,0);
  xpxi = (r'*r)\eye(nvar);
else % use Cholesky for very large problems
  xpxi = (x'*x)\eye(nvar);
end;

results.beta = xpxi*(x'*y);
results.yhat = x*results.beta;
results.resid = y - results.yhat;
sigu = results.resid'*results.resid;
results.sige = sigu/(nobs-nvar);
tmp = (results.sige)*(diag(xpxi));
sigb=sqrt(tmp);
results.bstd = sigb;
tcrit=-tdis_inv(.025,nobs);
results.bint=[results.beta-tcrit.*sigb, results.beta+tcrit.*sigb];
results.tstat = results.beta./(sqrt(tmp));
ym = y - mean(y);
rsqr1 = sigu;
rsqr2 = ym'*ym;
results.rsqr = 1.0 - rsqr1/rsqr2; % r-squared
rsqr1 = rsqr1/(nobs-nvar);
rsqr2 = rsqr2/(nobs-1.0);
if rsqr2 ~= 0
results.rbar = 1 - (rsqr1/rsqr2); % rbar-squared
else
    results.rbar = results.rsqr;
end;
ediff = results.resid(2:nobs) - results.resid(1:nobs-1);
results.dw = (ediff'*ediff)/sigu; % durbin-watson



function x = tdis_inv (p, a)
% PURPOSE: returns the inverse (quantile) at x of the t(n) distribution
%---------------------------------------------------
% USAGE: x = tdis_inv(p,n)
% where: p = a vector of probabilities 
%        n = a scalar dof parameter
%---------------------------------------------------
% RETURNS:
%        a vector of tinv at each element of x of the t(n) distribution      
% --------------------------------------------------
% SEE ALSO: tdis_cdf, tdis_rnd, tdis_pdf, tdis_prb
%---------------------------------------------------

%       Anders Holtsberg, 18-11-93
%       Copyright (c) Anders Holtsberg
   

  if (nargin ~= 2)
    error ('Wrong # of arguments to tdis_inv');
  end

s = p<0.5; 
p = p + (1-2*p).*s;
p = 1-(2*(1-p));
x = beta_inv(p,1/2,a/2);
x = x.*a./((1-x));
x = (1-2*s).*sqrt(x);

function x = beta_inv(p, a, b)
% PURPOSE: inverse of the cdf (quantile) of the beta(a,b) distribution
%--------------------------------------------------------------
% USAGE: x = beta_inv(p,a,b)
% where:   p = vector of probabilities
%          a = beta distribution parameter, a = scalar
%          b = beta distribution parameter  b = scalar
% NOTE: mean [beta(a,b)] = a/(a+b), variance = ab/((a+b)*(a+b)*(a+b+1))
%--------------------------------------------------------------
% RETURNS: x at each element of p for the beta(a,b) distribution
%--------------------------------------------------------------
% SEE ALSO: beta_d, beta_pdf, beta_inv, beta_rnd
%--------------------------------------------------------------

%       Anders Holtsberg, 18-11-93
%       Copyright (c) Anders Holtsberg
% documentation modified by LeSage to
% match the format of the econometrics toolbox

if (nargin ~= 3)
    error('Wrong # of arguments to beta_inv');
end
 
if any(any((a<=0)|(b<=0)))
   error('beta_inv parameter a or b is nonpositive');
end
if any(any(abs(2*p-1)>1))
   error('beta_inv: A probability should be 0<=p<=1');
end

x = a ./ (a+b);
dx = 1;
while any(any(abs(dx)>256*eps*max(x,1)))
   dx = (betainc(x,a,b) - p) ./ beta_pdf(x,a,b);
   x = x - dx;
   x = x + (dx - x) / 2 .* (x<0);
end
    

 
function pdf = beta_pdf(x, a, b)
% PURPOSE: pdf of the beta(a,b) distribution
%--------------------------------------------------------------
% USAGE: pdf = beta_pdf(x,a,b)
% where:   x = vector of components
%          a = beta distribution parameter, a = scalar
%          b = beta distribution parameter  b = scalar
% NOTE: mean[(beta(a,b)] = a/(a+b), variance = ab/((a+b)*(a+b)*(a+b+1))
%--------------------------------------------------------------
% RETURNS: pdf at each element of x of the beta(a,b) distribution
%--------------------------------------------------------------
% SEE ALSO: beta_d, beta_pdf, beta_inv, beta_rnd
%--------------------------------------------------------------

%       Anders Holtsberg, 18-11-93
%       Copyright (c) Anders Holtsberg
% documentation modified by LeSage to
% match the format of the econometrics toolbox
  

if (nargin ~=3)
    error('Wrong # of arguments to beta_pdf');
end

if any(any((a<=0)|(b<=0)))
   error('Parameter a or b is nonpositive');
end

I = find((x<0)|(x>1));

pdf = x.^(a-1) .* (1-x).^(b-1) ./ beta(a,b);
pdf(I) = 0*I;


function y = tdis_prb(x,n)
% PURPOSE: calculates t-probabilities for elements in x-vector 
%---------------------------------------------------
% USAGE: y = tdis_prb(x,n)
% where: x = vector containing computed t-values
%        n = degrees of freedom parameter
%---------------------------------------------------
% RETURNS:
%        y = a vector of marginal probability levels
% --------------------------------------------------
% SEE ALSO: fdis_prb(), chis_prb
%---------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jpl@jpl.econ.utoledo.edu


if nargin ~= 2; error('Wrong # of arguments to tdis_prb'); end;
if n <=0; error('dof is negative or zero in tdis_prb'); end;

x2 = n./(n+x.^2);
one = find(x2 >= 1);
if length(one) > 0
    x2(one,1) = 1-1e-12;
end;
zip = find(x2 <= 0);
if length(zip) > 0
    x2(zip,1) = 1e-12;
end;

tmp = 1.0 - 0.5*betainc(x2,0.5*n,0.5);
y = 2*(1-tmp);  


function xlag = mlag(x,n,init)
% PURPOSE: generates a matrix of n lags from a matrix (or vector)
%          containing a set of vectors (For use in var routines)
%---------------------------------------------------
% USAGE:     xlag = mlag(x,nlag)
%       or: xlag1 = mlag(x), which defaults to 1-lag
% where: x = a matrix (or vector), nobs x nvar
%     nlag = # of contiguous lags for each vector in x
%     init = (optional) scalar value to feed initial missing values
%            (default = 0)
%---------------------------------------------------
% RETURNS:
%        xlag = a matrix of lags (nobs x nvar*nlag)
%        x1(t-1), x1(t-2), ... x1(t-nlag), x2(t-1), ... x2(t-nlag) ...
% --------------------------------------------------
% SEE ALSO: lag() which works more conventionally
%---------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jpl@jpl.econ.utoledo.edu

if nargin ==1 
n = 1; % default value
init = 0;
elseif nargin == 2
init = 0;
end;

if nargin > 3
error('mlag: Wrong # of input arguments');
end;

[nobs, nvar] = size(x);

xlag = ones(nobs,nvar*n)*init;
icnt = 0;
for i=1:nvar;
for j=1:n;
xlag(j+1:nobs,icnt+j) = x(1:nobs-j,i);
end;
icnt = icnt+n;
end;

function z = trimr(x,n1,n2)
% PURPOSE: return a matrix (or vector) x stripped of the specified rows.
% -----------------------------------------------------
% USAGE: z = trimr(x,n1,n2)
% where: x = input matrix (or vector) (n x k)
%       n1 = first n1 rows to strip
%       n2 = last  n2 rows to strip
% NOTE: modeled after Gauss trimr function
% -----------------------------------------------------
% RETURNS: z = x(n1+1:n-n2,:)
% -----------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jpl@jpl.econ.utoledo.edu

  [n junk] = size(x);
  if (n1+n2) >= n; 
     error('Attempting to trim too much in trimr');
  end;
  h1 = n1+1;   
  h2 = n-n2;
  z = x(h1:h2,:);
  
